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Abstract. The classical method for deriving the macroscopic dynamics of a lattice 
Boltzmann system is to use a combination of different approximations and expan- 
sions. Usually a Chapman-Enskog analysis is performed, either on the continuous 
Boltzmann system, or its discrete velocity counterpart. Separately a discrete time 
approximation is introduced to the discrete velocity Boltzmann system, to achieve a 
practically useful approximation to the continuous system, for use in computation. 
Thereafter, with some additional arguments, the dynamics of the Chapman-Enskog 
expansion are linked to the discrete time system to produce the dynamics of the 
completely discrete scheme. 

In this paper we put forward a different route to the macroscopic dynamics. We 
begin with the system discrete in both velocity space and time. We hypothesize that 
the alternating steps of advection and relaxation, common to all lattice Boltzmann 
schemes, give rise to a slow invariant manifold. We perform a time step expansion 
of the discrete time dynamics using the invariance of the manifold. Finally we 
calculate the dynamics arising from this system. 

By choosing the fully discrete scheme as a starting point we avoid mixing ap- 
proximations and arrive at a general form of the microscopic dynamics up to the 
second order in the time step. We calculate the macroscopic dynamics of two com- 
monly used lattice schemes up to the first order, and hence find the precise form 
of the deviation from the Navier-Stokes equations in the dissipative term, arising 
from the discretization of velocity space. 

Finally we perform a short wave perturbation on the dynamics of these example 
systems, to find the necessary conditions for their stability. 

1 Introduction 

The Boltzmann Equation is a key tool within statistical mechanics, used to 
describe the time evolution of gases, and with some extensions, other flu- 
ids also. In this work we are concerned with calculating the dynamics of an 
ideal gas. This is achieved by calculating the statistical behaviour of sin- 
gle particles, that is the distribution of their positions in phase space. In 
particular by fixing a time point and integrating across velocity space, it 
is possible to calculate the macroscopic quantities of a fluid across space. 
Performing such an integration allows us to take a snapshot of the dynam- 
ics at any time point. If we are concerned, however, with discovering the 
rates at which the macroscopic variables change, we need to apply additional 
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techniques and assumptions to the Boltzmann Equation. A common choice 
of technique to derive these macroscopic dynamics is the Chapman-Enskog 
procedure |1|2) . This method involves calculating the dynamics of the distri- 
bution of the particles at different orders, within the Boltzmann Equation 
following a perturbation by a small parameter, the Knudsen number. Under 
such a perturbation the convective dynamics of the fluid will appear at the 
zero order, and the dissipative dynamics at the first order. The final result is 
that following such a treatment, the Navier-Stokes equations will be revealed. 
At the second-order and beyond additional terms give rise to Burnett and 
super-Burnet type equations repectively. 

Despite the Boltzmann Equation recovering the Navier-Stokes equations, 
to first order in the Knudsen number at least, a practical investigation into the 
macroscopic dynamics of the Boltzmann system is not neccesarily complete. 
In order to solve the Boltzmann equation numerically, a discretization of 
both time and velocity space is necessary. To be clear, in many cases for 
lattice-Boltzmann methods, a discretization of space is not necessary, as a 
good choice of velocity set and time step size can result in the Boltzmann 
equation being solved on a discrete lattice subgroup of space points, hence 
suffering no extra error in this regard. 

Within the discrete velocity, continuous time scheme the method of choice 
to evaluate the macroscopic dynamics of the system has remained to be the 
Chapman-Enskog procedure. This is presented in a number of different ways 
|3I4) . however the crux of the approach remains the expansion in the small 
parameter the Knudsen number. The necessary discretization of the velocity 
performed and described in some detail in Sec. El already gives rise to an 
additional error from the Navier-Stokes equations. Due to the approximation 
to the Maxwell distribution, the Navier-Stokes equations are only recovered 
exactly as the Mach number tends to zero. Where the Mach number is large, 
additional viscous (dissipative) terms appear [5] 

In order to move from continuous time to the lattice system, an Euler 
step is used to approximate the time derivative. Collisions should happen 
exactly once per time step, therefore the rate at which collisions occur is given 
by the time step itself. Furthermore, if the time step is small, in the same 
scale as the Knudsen number inherited from the continuous time system, the 
asymptotics of the order by order expansion may be compromised. Bearing 
in mind the complexity of combining the necessary number of terms from the 
Euler approximation, along with the existing expansions from the continuous 
system and taking into account the additional small parameter, we present, 
in this work, a possibly simpler route to the macroscopic dynamics. 

Our starting point is, in fact, not the continuous Boltzmann equation, 
but the discrete time system itself, in this sense we work in parallel to the 
historic lattice gas automata idea. We choose that the time step is small, and 
it is this parameter which we use for our asymptotic analysis. By choosing a 
discrete scheme such that the zero order dynamics give the Euler equations, 
we show in Section [3] that we retrieve the same computational system as the 
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discrete time Boltzmann anyway. We pursue the dispersive dynamics as the 
higher order dynamics, in the time step, of the difference scheme that we 
have chosen. Such a perspective is motivated by, for example Goodman and 
Lax [B] where it was shown that a particular difference scheme applied to the 
partial differential equation, 

d d 

—u(x,t)+u(x,t)—u(x,t) = 0, (1) 

recovers at the second order in the space difference parameter A, the KdV 
equation, 

d did 3 

—u{x,t)+u(x,t)—u(x,t) + -A 2 u(x,t)-^u(x,t) = (2) 

In tandem with the discrete time step asymptotics we hypothesize the ex- 
istence of a slow Invariant Manifold j7], we calculate the general form of this 
manifold and use it to find an expression for the macroscopic dynamics of 
general discrete time systems, with both discrete and continuous velocities. 
As we have stated our choice of the small parameter to be used in the asymp- 
totics is the discrete time step itself e. The dynamics of the quasi-equilibrium 
approximation (the Maxwell distribution in the continuum) define the zero 
order dynamics, higher order dynamics are given by the correction to the 
equilibrium of the same order [5] . We match the dynamics of the distribution 
function at microscopic and macroscopic [5] levels to find an expression for 
the first order non-equilibrium component of the distribution function. To- 
gether the zero and first order components of the distribution are sufficient 
to calculate the macroscopic dynamics up to the first (dissipative) order. 

As well as deriving a general expression of the macroscopic dynamics, we 
additionally provide two examples of both discrete and continuous velocity 
systems. We find that despite the qualitative difference between continuous 
and discrete time systems, in discrete time we can still recover the Navier- 
Stokes equations in a continuous velocity system. In the discrete velocity 
system, however, the discretizations examined display, as expected, additional 
errors in the dissipative part. The precise form of these errors is subject to 
the discrete velocity set chosen and given for the two common examples we 
use. 

Finally we test, by a short wave perturbation, the stability of the dynamics 
of the discrete velocity system. 



2 Background and Notation 

We begin with a short summary of the background to this work which will 
introduce some of the requisite ideas and notation. This brief discussion will 
be sufficient for this work, for many more details regarding the background 
and theory to Boltzmann systems a large number of works exist by several 
authors 03Enn|. 
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The Boltzmann Equation is concerned with the time evolution of what 
might be described physically as the density of particles. This density function 
is denoted / and is a function of space, velocity space and time f = f(x, v, t). 
It is given as follows; 

§- t f + vV x f = Q c (f)- (3) 

Already we need a certain amount of details to formalize what we have writ- 
ten. We begin with the space variable x, this is a smooth fc-dimensional 
manifold, for example Euclidean space, or a torus, in Unite dimensions. The 
velocity variable v ranges over the same space. In the continuous time system 
above t is a single real variable, however we will dispense with this very soon. 
The final notation to mention in [3] is the collision integral Q c , a differentiable 
transformation of the population function. More detailed properties which 
we require of this function will be given later. 

In fact we are not at all concerned with the continuous time Boltzmann 
evolution, but the corresponding discrete time system: 

f(x + ev, v,t + e) = f(x, t) + Q(f(x, t)). (4) 

This is the discrete time Boltzmann system which is at the core of this work. 
In contrast with the continuous time system we have introduced the time 
step e, which will play a key role in our analysis. The time step should be 
small and it restricts admissible values of time to the subgroup eZ € R. 

The left two terms in Eq.@]are collectively termed advection or free flight, 
if the collision integral is omitted, the exact evolution of the population func- 
tion at a point may be written as f(x, v,t + e) = f(x — ev, v, t). The physical 
interpretation is that particles move freely under their own momentum with 
no interaction between themselves. For our analysis we can use the fact that 
advection is smooth and we will be able to take a Taylor series expansion to 
any finite order that we need. The rightmost term is again called the colli- 
sion integral and is denoted slightly differently to notate that it may differ 
somewhat in form from its continuous counterpart. 

Rewriting the right hand side of Eq. |4] leads us to a perspective which 
will be key to our approach. Rewriting the collision part as F(f) = f + Q(f), 
allows us to present Eq. @] differently: 

f(x,v,t + e)=F(f(x-ev,v,t)). (5) 

In this way, the discrete time Boltzmann system can then be thought of 
as a superposition of the advection and collision operations. Advection and 
collision steps happen in turn, this is qualitatively different from the con- 
tinuous time system where both operations might be said to be happening 
simultaneously at all times. 

The evolution of the particle density function / can be said to describe 
the microscopic evolution of the system. To recover the macroscopic system 
we need to sum accross the velocity space. Here then is the key distinction 
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between the continuous (v e R fc ) and discrete (v <E {v a }a = 0...n) veloc- 
ity systems. In the discrete velocity case, instead of considering the density 
function to be a function of velocity, we define different density functions /" , 
one per velocity vector and we can denote then by / the ordered set of them 
/ = f a (x, t), a = 0, . . . , n. For example the density of the fluid is calculated 
in the continuous and discrete schemes respectively as 



In fact which macroscopic variables M we use are not important to our initial 
analysis. The properties that are important for the macroscopic moments are 
firstly that the operator to, which recovers the macroscopic moments from the 
density function / (to(/) = M) is linear. The second property is that these 
moments are invariant under the collision operation, that is m(/) = m(F(f)). 

Defining the macroscopic variables allows us to discuss another property 
of the collision operation F. The quasi-equilibrium is the unique vector / = 
such that F(f^) = Jm and m(flf) — M. Existence and uniqueness of 
this quasi-equilibrium will be assumed in this analysis. There is exactly one 
equilibrium point per value of M and together they form the quasi-equilbrium 
manifold through the space of the density function. 

3 The discretization of Velocity Space 

While we calculate the dynamics of the continuous velocity system later, they 
serve only for the purpose of calculating the error incurred by an approxi- 
mation. For practical computations we would use a discrete velocity system. 
What we have control of is which discrete approximation to use. Two different 
approaches will lead us to the same result. 

Firstly we consider an actual discretization of the Maxwell-Boltzmann dis- 
tribution in D dimensions, that is the known quasi-equilibrium in continuous 
velocity space. 



The 2 + D thermodynamic moments here are the density p, momentum u and 
temperature T . This distribution is multiplied with low order polynomials 
of the velocity and integrated to retrieve the macroscopic moments. 



(6) 






(8) 



We can view the velocity set as a quadrature approximation to these integrals 
[?] . With such a perspective the choice of which quadrature to use may seem 
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obvious, integrals of Gaussian type (as are the moment integrals) can be 
integrated exactly using a Gauss-Hermite quadrature. Unfortunately several 
stumbling blocks prevent us from performing this exact integration. 

The first problem is the necessary change of variable in the equilibrium 
distribution. For a Hermitian quadrature of any order, the nodes should be 
distributed symmetrically about the centre of the Gaussian (in this case u), 
and the coordinate of integration should be normalized. Effectively it is nec- 
essary to apply a change of variable to the moment integrals so that the 
exponential term is of the simple form exp(— v 2 ). Unfortunately, to do this 
in practice would require us to know both u and T before performing the 
integration. Of course we may still be able to solve such a system, perhaps 
by some iterative method (moving the quadrature nodes) up to an arbitrary 
degree of accuracy, but only if we can evaluate the density function anywhere 
we choose. In practice of course we can only store a finite number of eval- 
uations of the density function, and these points of evaluation are chosen 
pre-emptively to be the same across all lattice sites and time steps (in order 
that advection is an exact operation) . Because of these factors we should not 
expect to be able to integrate exactly a density function of Gaussian type 
with a general momentum and temperature. 

The popular method to partially rectify this problem is to assume that u 
is close to zero. If we choose to work in an athermal system (where T is some 
constant) then we can exand the Maxwell distribution about the point u = 
up to the second order giving us 

f-v 2 \ ( vu (T-v 2 )u 2 \ 
r = ,(2flT)- exp ^— J (l + - - [ —^) 0) 

The second order expansion is taken to get sufficient terms that the temper- 
ature moment integral is calculated correctly. With a constant T , and u no 
longer affecting the midpoint of the distribution, the same nodes and weights 
can be used for all space points and time steps. Assuming the proper trans- 
formation of variable these nodes and weights are the standard hermitian 
ones. 

The alternative method of defining the quasi-equilibrium is by solution of 
a linear system. Since the moments are calculated by linear combinations of 
the discrete quasiequilibria, these equilibria can be found by the inversion of 
a matrix of components of the discrete velocities, to illustrate this we give an 
extremely simple example. We consider an athermal one dimensional system 
with three discrete velocities {—1, 0, 1} and fixed temperature 1/3. We create 
a matrix where each line represents an elementwise power of the velocities 
(up to the second order) and the corresponding vector of moments we desire. 
Together this creates the following linear system. 



"ill" 




p 


-1 1 


• r = 


pu 


1 01 




_p/3 + pu 2 _ 
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Solution of the above system yields a vector of equilibria 

/ cq = 1 1 (P - ipu + 3pu 2 ) ,^(p- ^-pu 2 ^ , I (p + 3pu + 3pu 2 ) | (11) 

Inspection reveals this to be exactly the same, with weights pre-included, as 
the discretized Taylor approximation to the Boltzmann (Eq. [9]) where the 
same velocity set is applied. We note that in this example we have exactly 
the necessary amount of velocities (that is one per moment) and hence have a 
square matrix. As in the case of any linear system, were we to have fewer ve- 
locities it would become likely that the system would become unsolvable. For 
each additional discrete velocity that we add above the necessary amount, we 
can impose an additional linear constraint. Often the choice would be to zero 
higher order moments of the equilibrium distribution, that is enforcing that 
the sum of the equilibria multiplied with polynomials of the discrete velocity 
components, with order higher than two, should be zero. Another option is 
to create additional, non-hydrodynamic moments in order to suppress insta- 
bilities [12]. 

4 Invariant Manifolds for Discrete Time Boltzmann 
Systems 

We will, by finding a general form for an invariant manifold, calculate general 
microscopic dynamics for discrete time Boltzmann systems. To do this we 
need to make some assumptions regarding the stability of collisions and the 
smoothness of the distribution function. 

For the next subsections we introduce the notation that fields of distri- 
bution functions of moments will be written in gothic font, hence the field of 
macroscopic variables, for example, will be denoted 

We assume that collisions are stable, and for any admissible initial state / 
iterations F p (f) converge to unique equilibrium point / cq exponentially fast 
and uniformly: 

ll^(/) " C(/)H < Cexp(-Xp)\\f C (f) \\, (12) 

where the Lyapunov exponent A > and pre-factor C > are the same for all 
admissible /. In the limit e = there is no free flight, the field of macroscopic 
variables DJt does not change, and the field of distributions f converges to the 
local equilibrium field by repeated application of the collision operation. 
Since each collision occurs instantaneously, the superposed collisions become 
a projection onto the local equilibrium in zero time. 

In order to discuss small e > we need to evaluate the change of macro- 
scopic variables in free flight during time step e. To find the qth term of the 
non-equilibrium density function of a discrete velocity system, we make two 
assumptions: 
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— The first q derivatives of /m along vector fields v a exist and are bounded. 

— The differentials (-D«/m)i are uniformly bounded for all M. 

Our expression for the manifold will be of the form of an aysmptotic 
expansion in the small parameter e, 

/ inv = / (0) +e/ (l) +e/ (2) +o(e) _ (13) 

Our first goal is to find a prescription for this term. The zero order term of 
this is simply given by the quasi-equilibrium distributions, that is = / cq . 
For the first order term we will take expansions of the distribution function 
in terms of time and of moments and equate them. That is for each order 
in epsilon we can take the effect of a complete LBM step (advection and 
collision) and match the effect on the distribution function to that of taking 
the Taylor approximation of the manifold through the moments up to the 
same order. In other words we match the dynamics of the microscopic and 
macroscopic scales on an order by order basis. 

4.1 The Invariance Equation 

The procedure we use can also be described in terms of the invariant manifold 
hypothesis. Coupled steps of advection and collision form a chain of states 
of the population function belonging to the manifold. Since the number of 
discrete velocities used is normally larger than the number of macroscopic 
moments, there are an infinite number of possible population distributions 
which can give rise to the same configuration of moments, however only one of 
these distributions exists on the manifold. We use a Taylor approximation to 
the manifold and match it with a single coupled step to find the components, 
at different orders of the time step, of the distribution function. 

If we consider f 1 ^ to be the field of population distributions on the mani- 
fold with corresponding field of macroscopic moments 9Jt then in a continuous 
time system this invariance property can be defined as 

D t f^ v = D m f™ ■ DM- (14) 

Here the derivative Dot indicates the derivative through the field of macro- 
scopic moments 9Jt, which the field of distributions on the manifold fj^ v is 
parameterized by. Altogether the rate of change of the population function is 
equal to the rate of change of the moments multiplied by the change of the 
populations with respect to the moments. The discrete time analogy of this 
is given by, 

(fim)' = & (15) 
where the prime notates the next time step, therefore the left hand side of 
this equation is given by Eq. 

From now on we dispense with the gothic notation for fields of variables 
and simply use a standard font to describe distribution functions or macro- 
scopic variables. 
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4.2 The expansion of the distribution function following a step 
in the LBM chain 

A little abuse of notation will make the same calculations in this section both 
brief and meaningful. No distinction will be made here between continuous 
and discrete velocity systems, however implicitly there is one. As previously 
/ = f(x, v, t) for the continuous case, and / = f a (x,t), a = 0, . . . , n in the 
discrete case. For the purposes of our analysis, in both the continuous and 
discrete velocity systems the space coordinate x is continuous. Due to this, 
even though we may choose to only evaluate / at a discrete number of lattice 
sites, the space derivative D x f remains well defined, even for the discrete 
velocity case. The notation is similar for moment derivatives, where Dm/m 
refers to the change in moments of the function / evaluated at the moments 
given by the subscript to /. When the moment operator m is applied this 
refers to the integral in the continuous case and the summation in the discrete 
case as given in Eq. [6] 

With these ideas in mind the procedure below is valid for both cases. 
The first ingredient for the time step expansion is the Taylor series of the 
advection operation up to the required order in e. For the first order we have 

f(x-ev) = f-evD x f + o(e). (16) 

Combining this with ()13[) we have to the first order, 

f{x - ev) = /(°) - ev ■ D x /(°) + e/« + o(e). (17) 

Applying a collision operation gives the complete, composite discrete time 
step, 

U'm)' = F (/ (0) - ev ■ D B /<°> + ef W + o(e) ) . (i 8 ) 

The second ingredient is to use a linearised version of the collision operation, 
this is sufficient to get the first order populations correctly. Here the lineari- 
sation is made about the equilibrium corresponding to the populations to be 
collided, 

Due to the linearity we can move the error term in Eq [18] outside the collision 
altogether. The linearisation is then made about the equilibrium defined by 
the moments of the first order advected populations 

M[ = m (/<°> - ev ■ D x f^ = M + m (-«; • D x ff^ , (20) 

Finally then for the first order approximation to the next step through the 
time step expansion we have, 

(Pm)' = ft. + (D f F) n> (/$ + e/# - ev ■ D w fjg> - /£,) + o(e). (21) 
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4.3 The expansion of the invariance equation following a time 
step 

With the expansion of the left hand of (fT5")) complete we consider the right 
hand side. Here we find the Taylor expansion of the invariant manifold up to 
the linear term so, 

fw = !m + (DmJm) ■ m {-ev ■ D x f M ) + o(e). (22) 
Substituting (fT3"|) into (|2"2"]| we have 

/m< - C + ef$ + (DmC) ■ m (-ev ■ D x f$) + o(e). (23) 
We can now equate (|2"Tj) and (j2"3")l for a first order approximation to (|T5|) . 

= /m + e/ff + (Atf/J?) • m (-ev ■ D x fj») . (24) 

Of course in a similar style to (f2"2"j) , 

= /i? + • m (-a, • D x /<?) . (25) 

Substituting back into (|2"4")l we have 

(DfF) ^.cq (e/« - eu • D./J? - (D M f$) ■ m (-ev ■ D x fj>?)) = e/« 

(26) 

This equation forms the prototype to find f$ for different possible collision 
operations, it implicitly gives the first order approximation to the invariance 
equation (fT5)l . It depends on the choice of the velocity set, the quasiequilib- 
rium and the collision integral. 

4.4 Example First Order Invariant Manifolds 

We consider two possible examples of collisions. The first example is the 
simple Ehrenfest step [T3] . 

F(f) = C u) (27) 

we immediately have, 

(D f F)f m(f) (f) = 0. (28) 

Substituting back into (|2"6")l , 

= 0- (29) 

This of course expected since using Ehrenfests steps for the collisions we 
should expect to return at every time step to the quasi-equilibrium manifold. 



The Invariant Manifold Approach to LBMs 11 
The second example of a collision operator is the BGK collision |14) . 

*"(/) = / + <"(C(/)-/)- ( 3 °) 

Differentiating we have 

(D f F) n (/) = (l-c)-/ (31) 
Substituting this into (|2"6jl . 
(1 - a;) - (e/« ev ■ D x f$ - (D M f$) ■ m (-ev ■ D x f$)) = e /£>. (32) 

We can multiply out (l32t and solve for, /J£ . 

T^f$ = {-v • D x f^) - (^M/|?) ■ m (-v ■ D x f$) (33) 

Figure Q] graphically demonstrates the f$ for the BGK collision type. In 
particular for this example there is a critical parameter value at w = 1. For 
w = 1 we recover the Ehrenfest step, for w > 1 we have the normal BGK 
over-relaxation where both of the coupled steps of advection and collision 
cross the quasiequilibrium manifold, one in each direction. 



4.5 Second Order Manifolds and an Example 

The next goal is to find an equation analagous to (|26|) for the second order 
term of the invariant manifold. During the next section we use a linear colli- 
sion operation, in this case the linearised collision we use produces the exact 
same result as the original collision. We restart the procedure using second 
order expansions where appropriate, the first of these is the Taylor expansion 
of the advected populations. 

fix -ev) = f-ev D x f + ^yD x (v DJ) + o(e 2 ) (34) 
The second order population expansion is also used. 

/in v = /( 0) + £/ (l) + e 2 /( 2) + o(e 2 } _ (35) 

Altogether the second order expansion of the advected populations is, 



f{x - ev) = / (0) - ev ■ D x f^ + *-v ■ D x (y ■ D x f®) 

+ e/ (i)_ e VAj(i) + e 2 /( 2) +o(e2) (36) 
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Fig. 1. Graphical representation of the for the BGK collision. Adding 

the fff term to the quasiequilibrium manifold gives the invariant manifold to 
first order in e. In particular the collision parameter u> is critical, for u > 1 the 
direction of the f$ term is inverted and consequently the invariant manifold 
is below (in the sense of this illustration) the quasiequilibrium. Therefore 
at each step the advection operation crosses the quasiequilibrium and the 
collision returns below it. 

We use the linearized collision integral in replacement of the original collision 
operation, 

(fu)' = fll + ( D f F )Q, (/m - ™ ■ D x f { V + € lvD x (v- D x f$) 

+ ef$ - e\ ■ D x f$ + e 2 /^ - /£,) + o( e 2 ). (37) 

where M 2 are the moments of the post advection populations to second order, 

M' 2 = m (f$ - ev ■ D x f$ + e l v -D x {v- D x f$) ~ e 2 v ■ D x f^ , (38) 

For the right hand side of (p~5]) we use a second order approximation to the 
invariant manifold, 

f$ = f^{^v a ) + AM 2 -DMf^ + \AM 2 -D M {AM 2 -D M fTn + o((AM 2 ) 2 ). 

(39) 
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where for berevity AMi is the difference between the moments of the post 
and pre advection populations to second order, 



AM 2 = m (-ev ■ D x f$ + ^v-D x (v D x f { °^ - e 2 v ■ D x f^j = M' 2 -M. 

(40) 

Substituting ((35|) and (gO) into ([39l we have 



& = $ + <f2? + <?f$ + m( ev • D m $ + ^vD x (v D x f<$) 
- e*v ■ D x f$) ■ D„fj® 
+ \m (-ev ■ D x f$) ■ D M (m (-ev ■ D x f$} ■ £> M /|?) 
+ em (-ev ■ As/]?) ■ D M f$ + o{e 2 ). 

(41) 

With the expansions of both sides complete we can equate (|37|) and (|41~|) . 



(i) 

M 



III + ( D f F )Q, [fu - ™ • D*f® + jv-D x (v D x f$) + ef 



- e 2 v ■ D f (1) + e 2 f {2) - f cq 

AO), ,(1), 2,(2) 
J M + t J M + € JM 

+ m (-ev ■ D x f$ + € -v-D x (v ■ D x f$) - e 2 v ■ D x f M } ) ■ D M f$ 



+ l m (- ev . D x f$) ■ D M (m (-ev ■ D x f$) ■ D M f 
+ em (-ev ■ D x f^ ■ D M f 



;(1) 
' M 



(42) 



Analagously to the first order case we note that, 

fl\ =/$ + m( ev ■ D x fj» + e lvD x (v D B /<°>) - 

e 2 v-D x f^)-D M f^ (43) 
+ m (-ev ■ D x f$) ■ D M (m (-ev ■ D x f^) • D M f$) ■ 
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Substituting this back into (fl2l we have the final prototype for f M ' which 
this time is given implicitly, 

{ d S f )q, (-ev-D x f$ + ^v-D x (v ■ D x f$)+ef£> -e\-D x f$ +e 2 f$ 

- m (-ev ■ D x f$ + 6 -vD x (y- D x f$) - e 2 v ■ D x f$^j ■ D M f$ 

- \m (-ev ■ D s /<°>) ■ D M (m (-ev ■ D x ff?) ■ D M f { ^) ^ 

= *f$ + t 2 fM + em (-ev ■ D x f$) ■ D M f$ (44) 

We return to the BGK collision for a specific example of an f M term. (|44p 
becomes, 




Rearranging and equating terms with e order 2 gives us, 




We will not use these populations in the examples of macroscopic dynamics 
which we calculate in the next section. For the dissipative dynamics the 
first order populations are sufficient. We expect that this second order part 
should give rise to macroscopic dynamics equating to the Burnett equations 
in a continuous velocity system, with some additional error terms if a discrete 
velocity set is used. 
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5 Macroscopic Equations 

In this section we are concerned with deriving equations for the macroscopic 
dynamics arising from several different example lattices. We expect that the 
lattice parameter e should partly govern these dynamics and that the 1st 
order macroscopic dynamics should be governed by the 1st order population 
functions. 

In order to find these dynamics we project the microscopic flow (advec- 
tion) up to the required order, following one time step, onto the invariant 
manifold up to the same order [7|15j . 

We can immediately perform a Taylor expansion in time on the macro- 
scopic dynamics, 

M' = M + e^-M + o(e) (47) 
at 

We expect that the final model should be given in terms of a time derivative 
of the moments, we write this in a power series in terms of e, 

— M = V {Q) + eW (1) + o(e) (48) 

Combining these two we have 

M' = M + e<f (0) + o(e) (49) 
Equating d^HJ) and (gH]) we have 

* (0) = m (-v ■ D w f$) (50) 
The corresponding second order approximation of the moments in time is 

M' = M + e(^ (0) + e<f (1) ) + € - J^°) + o(e 2 ) (51) 
Equating terms on the second order of e we have, 

m Q« • D x (y ■ D x f$fj +m(-v D x f$) = ^\t) + ~* i0) (t) (52) 
or 

* (1) (*) = m (lv ■ D x (y ■ D x f$)~) + m (-v ■ D x f$) - \^ ( °\t) (53) 

Although in this notation the evolution of the macroscopic moments of the 
continuous and discrete velocity systems can be written in the same way, in 
practical examples we have no reason to suspect that they should be necessar- 
ily equivalent. In the next two sections we calculate some example, equivalent, 
discrete and continuous systems, in order to compare them. 
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6 Discrete Velocity Examples 

We will now demonstrate the exact first order dynamics of a popular choice of 
lattice scheme in one and two dimensions. The athermal schemes we consider 
are typically described in shorthand by the dimension within which they 
operate and the number of velocities used to form the lattice in the form 
DnQm, where m and n are integers representing the number of dimensions 
and velocities respectively. The general quasi-equilibrium for these systems, 
including the two examples we use can be written in a general form 

f a < (0) =W p(l + — + ' U? —) (54) 

This equilibrium defines an athermal system where the temperature is fixed. 
To complete the definition of the discrete system requires only the selection 
of a velocity set and some accompanying weights W a . 

6.1 An athermal three velocity lattice (D1Q3) 

Our 1-D example lattice is one of the most common, the athermal 1-D lattice 
with 3 velocities. In this example the velocity vectors are {—1,0, 1} and the 
speed of sound c s = l/y/3 hence the equilbrium populations are derived from 
the general formula for athermal quasi-equilibria in any dimension where the 
additional parameters the weights W a are { g , | , ^ } . 
For this case the populations are, 



1 



- jp(l -3u + 3< 2p \1 - — , p(l + 3u + 3u 2 ) | . (55) 

For this lattice with unit distances we note that v 4 = v 2 , v 1 — v 3 etc. We 
calculate the two components of using the formulas for the moments. 
We have for the density derivative, 



and for the momentum derivative 

!_0_ f a,(O) 

dx jM dx 



4 0) = -E^(- a ) 2 |:/M (0) = -^E w ^ a ) 2 fM {0) 



Now we examine the individual moments of the first order part in the case 
of the one dimensional lattice, as before we begin with the density, 
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The second term here is the space derivative of the momentum of the fff 
which equals zero due to all moments of non equilibrium components be- 
ing zero and the first term can be calculated immediately from the quasi- 
cquilibrium. Therefore, 

The time derivative of tf^ can be calculated by the chain rule, 



dt 1 d P dt + d P u dt dx 2 dxAs +p J { ' 



Substituting this back in we have, 



\&_(p_ (puf \ 1 d 2 (p 
~ 2dx 2 \3 + p J 2dx 2 



^+pu 2 )=0. (61) 



For the momentum moment we have 

d 2 „„.foi v^, 9 „„. m 1 9 



V. 



(i) 



2 



Recalling that that v 3 = v 1 we can simplify the first term so, 

1 W a\3 ® ta,(0) 1 Y^ a ra,(Q) 1 9 2 

- 2 U V ) = 2 J> - 2^"' (63) 

a a 

For the second term we need to calculate the terms. To do this we need 

to specify a collision type, we use the BGK collision described above (|33|) . 

w m 9/9 , 2 \ 9 9 2 

9 A 3 2 \ 9 n 9 2 

2w— p + (-2 - 3u 2 ) — pu + 6u— pu 2 ] . (64) 
ox ox ox I 



This gives then 



v^/ ax2 d Wi) I — a; 9 /2 9 / 2 2 \ 9 9 2 \ 



(65) 



Again using the chain rule, 



0t J 9p 9* dpu dt dx \ 3 / 1 9x 

9/2 9 /l 2 \ 9 9 , . 

^UW+U - " Jfe' u + at W u )■ (66) 
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Substituting this all back in we have, 



2 d (2 d ( 2 2 \ n d . 
—u—p + — — — u -z—pu + 2u—pu' 



uj-1 1\ d (2 d ( 2 ,\ d d 



2u dx \ 3 dx \ 3 / dx dx 

-2d 











-p(^ 2 -l) 






dx ) 



2uj dx 

The moment gradients are then to first order in e, 

d_ _ d_ 

dt P dx P 

d 0/1 2 \ 2-lo d ( 3 d / 2 2\ d , 

0^= U p+/Mr J ~ e ^r^ r fe p+p l 3u " 3 d- x u) 



(67) 



(68) 



6.2 An athermal nine velocity model (D2Q9) 

The 2-D example we consider is a popular 2d lattice consisting of 9 different 
velocities. If we identify v\ as the horizontal component of a vector and vi 
the vertical component then the set of velocities is 

o) ' (o) ' (l ) ' ( ) ' (-1 



(69) 



The equilibrium is then given by the polynomial formula (|54[) with corre- 
sponding weights 

(411111 1 1 11 

a ~ \9' 9' 9' 9' 9' 36' 36' 36' 36 J ( } 

As before we calculate the components of using the formulas for the 
moments although this time we have two momentum density momentums for 
the two dimensions. We have for the density derivative, 



- dxi 2^ v ^M q 2-~i 2 * M ( 71 ) 

a 

d d 
dx\ dx2 
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for the first momentum derivative 

a 

--^Z.( w i)/m ~^E^/m (72) 

a a 

9/1 2 \ a 



and for the second momentum derivative 

^ 0) =E^(-^-^/m 0) ) 

a 

--t!:E«/S (0) -^-E^) 2 /S (0) (73) 



9 9/1 



= --K—puiu 2 - ^— -p + pu 2 

OX\ 0X2 \3 

The first order density moment is given by, 

= \ E • *>* (« ■ ^)) + E (-« • - ^ (0) - ( 74 ) 

Again we observe that the second term is the space gradient multiplied with 
the momentum densities of the first order populations and hence is zero, for 
the first term we have 



J2(vD x (v- D x f 



V ( — (v a ) 2 f aM 4- 2 ^ ,.a,.a f a.(0) , J?l(.,aN2 f a,(0) 

9 2 /l 2 \ _ d 2 d 2 /l 



= dx\ U' + J + 2 dx^dx-/ UlU2 + dxJ^ P + pU2 >' (?5) 



and for the third term 



dt 1 dp dt dpui dt dpu 2 dt dx\ 2 dx 2 3 

0/0/1 2 \ 

-^(-^^-^G p+pu *ii- (76) 

hence subtracting these we have = 0. 
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For the first second order momentum density we have 



4 1] = \ £«? (v- D x {v. U./J°>)) + £«f • D./£>) - i^f. 

(77) 

Examining each term in turn more closely we have for the first term 



a 



(0) 







9 



9x 




1 .9 



1 







a — I Pa — "i + "i" 

oxi \ ax\ o^i o «j / 

d (i d id id l d \ 

+ ?P^— U2 + q"2t^P + u i + n"!^-/ 1 ■ (78) 



The first order populations are given in Appendix [5J these give us for the 
second term 



(i) 

u 2 J 

2 



£J - 1 / d ( o d 2 9 / 2 2\d 

v aa;i \ rai 0x2 \ o / crai 

9 2 d 

+2pUlU 2 - Ml + pM x - M 2 

aa; 2 aa;2 

d ( 2 9 2 ^ r, <3 r> ^ 
"5 M^M2 "7: p + MlM 2 7^ P + 2 j OMiM27^ Ml + 2/9MiM 2 - M 2 

0x2 \ 0x1 ox 2 ax\ 0x2 

2 1 \ d l / , 1 \ d 



pu -~r dx- 2 ui + { pUi ~r dx~, U2 



(79) 
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and for the third term 



a^o) = cw^dp + m^dpu, + d^ v> d P u 2 

dt 2 dp dt dpui dt dpu 2 dt 

+ (- J sr»-w-)<-^ 
= {( ui + d^ p + [r 2 + u ^ 2 ) d^ p 

„ 2 \ d fl 2 \ d n d \ 

+ U p + pUi ) d^i ui + U p + pUi ) d^ U2 + 2puiU2 o^ ui ) 

d ((\ 2 \ d fi 2 \ d 

+ o^ 2 { [r 2 + u ^ 2 ) + U Ul + UlU2 ) d^ p 

d d d d \ 

+2pmu 2 - — mi + pu\- — mi + pu\- — u 2 + 2pmu 2 - — M 2 . 

OXi OX 2 OX\ ox 2 J 

(80) 



Combining all three terms we have, 



T m iuj-l 1\ / d ( , d 2 d / o 2 \ d 

9 2 5 
+2/3M1M2T; Ml + pu x - M 2 

did d d d 

+ — u\u 2 - — p + uiuj- — p + 2puiu 2 - Ml + 2pu\u 2 - — u 2 

dx 2 \ ox\ ox 2 axi ox 2 

1 \ d ( o 1 \ d 

' -)) 

(81) 



+ l^-oPl^« 1 +lp Ul --pl ^M 2 
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The final macroscopic equations for this particular lattice and quasiequilib- 
rium then are to first order 

d n 
—p = -D x pu 

d d fp 2 \ d 



dt^ 1 dx, V3 1 r L J dx 



C2 

2-cj / 8 ( n d 2 d / o 2 \ d 

d 2 d 

+2pu 1 u 2 -^—ui + pu 1 -^-u 2 

did d d d 

+ q — u\u 2 - — p + u x u\- — p + 2pmU 2 - — Ul + 2pmu 2 - — u 2 
ox 2 \ ax\ ox 2 axi ox 2 

2 1 \ d ( 2 1 \ d 

pu ^r) dx- 2 ui+ { pu ^3 p ) dx-^ 2 

(82) 

The second momentum density is available easily through symmetry. 



7 Continuous Velocity Examples 

We are now concerned with calculating, via the invariant manifold popula- 
tions, the macroscopic moments approximated by the LBM chain in a contin- 
uous velocity system. We select two examples, chosen to match the previous 
discrete velocity schemes. 



7.1 The athermal 1-D model 

The first continuous velocity model we will examine is one chosen to match 
the zero order dynamics of the discrete model studied in section [6.11 the one 
dimensional system with the three discrete velocities { — 1,0,1}. The contin- 
uous population function acting as the quasi-equilibrium is a specific case of 
the Maxwell distribution where the temperature is fixed, in this case to 1/3. 

/(o)=p ^ exp H (w - u)2 ) (83) 

With such a system the macroscopic variables are calculated as integrals 
rather than the sums in the discrete case. 

oo 

f<®dV = P 

-oo 

o 

vf {0) dv = pu (84) 

K) 

v 2 f {Q) dv = \p + pu 2 
3 
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Clearly this matches the moments retrieved in the discrete velocity case. Due 
to this we can, analagously to the discrete case, immediately write down the 
zero order macroscopic dynamics following equation 1501 

(86) 

In order to calculate the first order moments we expect that we shall require 
the first order continuous populations. These are also derived exactly as in 
the discrete case with the replacement of the sum, by the integral, in the 
calculation of the moments. Since we replicate the discrete case the collision 
we select is again the BGK collision and we derive the first order populations 
from equation l33l 

' = P\/?exp (-hv - u)A ■ (1 - 3« 2 + Qvu - 3 U 2 ) • ( #u 



1 - u r V 2vr V 2 V / \dx 

(87) 

Again we calculate the first order moments from the template given by equa- 
tion [53] 



n = iL ' [a^ 1 ) dv -L 'lsH*"i*' 

Exactly as the discrete case the second term here is the space derivative of the 
momentum of the f^ 1 ' which equals zero due to all moments of non equilib- 
rium components being zero and the first term can be calculated immediately 
from the quasi-equilibrium therefore, 

Again the time derivative of can be calculated by the chain rule, 

dt 1 ~ dp dt + dpu dt -~dx^ dxA^ PU ) m 
Substituting we have, 



2dx 2 V3 r J 2dx 2 V3 
For the continuous velocity momentum moment we have 
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Rearranging and performing the first two integrals gives us 

f'-w^^-K-IKs-))-^ <»> 

Again using the chain rule, this term is exactly as in the discrete case, 



1*1°) _d^d_p dW^d_pu_ _d_(l_ \ (0) _ d_ (o) 



d [2 d fl 2 \ d n d o . 

= ^U u ^ p+ U- u J^ +2u ^ ' (94) 

Substituting this all back in we have, 



2 V w 2yfeV 3 p fe y 2w dx \ rdx 

The moment gradients are then, for the continuous velocity system, to first 
order in e, 

d d 

d d fl ,\ 2-uj d / 2 d \ . . [yo) 

m pu= -Tx {r +pu ) - e —Tx \-rTx u ) +o(e) 

We immediately observe that several of the dissipative terms that appeared 
in the discrete velocity system do not occur when we use continuous velocities 



7.2 The athermal 2D model 

The next continuous velocity model we examine is the widely used athermal 
2d model. Again we use a specific choice of the Maxwellian distribution which 
matches the zero order moments given by the discrete velocity set. 

/(0) = p i cxp ("I ( {vi ~ ui)2 + {V2 ~ U2)2) ) (97) 

Again macroscopic variables are calculated by integrals over velocity space 



Jm 2 



/ f^dv 
Jm 2 


= p 


[ Vl f {a) dv 


= pui 


Jm 2 




[ v 2 f {0) dv 
Jm. 2 


= pu 2 


(vl+vDf^dv 


2 



(98) 
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Again we calculate the zero order moments, 



^ (0) = [ -v D x r M (0) dv 



' } f I v 2 rj° ] dv (99) 



dx, ./:/ •' " 8x 2 I 

8 8 



= — pui - 3 — pu 2 
ax\ ox 2 

and for the first momentum derivative 



vi (-v ■ Ac/S (0) J dv 
j£- I vlftfUv v 1 v 2 ftf ) dv (100) 



8(1 2 \ d 

for the second momentum derivative 
,(0) _ 



/ v 2 (-vD x f^°Adv 



f vir^dv (101) 



n 



dxi J R 2 M dx 2 

8 8/1 2 

= __ pUlW2 __^_ p + / , U2 

We again calculate the first order populations following equation! 
-/ ( 



w f (i) 



1 - U 

p ~h exp (~I ~ u ^ 2 + ^ 2 ~ u2 - )2 ^) 

1 — 3u? + 6«iui — 3u?) — — Mi 

OTi 

+ (— 3viV2 + 3viu 2 + 3v 2 ui - 3itiit 2 ) — ui 

8x2 

8 

+ (— 3viV2 + 3viu 2 + 3v 2 ui - 3ititt 2 ) — "2 

ax\ 

+ (l - 3uo + 6v 2 u 2 — 3u?,) — 112 
x ' 8x2 



The first order density moment is given by, 



(102) 



*™ = \ • Dx ( v ■ Dxf{0) ) dv ~L v ' Dxj ^ dv - 1^ 0) f m s 
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Performing the integrals of the first two terms we note that the second term 
is again zero therefore 



m 1 2 /l 2 \ 2 1 d 2 (\ 2 \ 10 r(0) 

(104) 

and exactly as in the discrete velocity system we have for the third term 



5^(0 = d^_ dp + cw^dpui + d^dpu2 = _ _ A^(Q) 

dt 1 dp dt dpu\ dt dpu2 dt dx\ 2 dx 2 3 



0/0/1 2 \ 



0a;i \ dxi \3^^^ 1 ) dx2 P 1 2 
0/0 0/1 2 

-0^l~0^ m ~0^U^ + ^ 2 



hence subtracting these we have =0. 



For the first second order momentum density we have 



(105) 



& = \ I * („ .d.( v . d./<«)) dv+ 1 Vl (_ • Dmf <») dv -\^ 

(106) 

Again performing the integrations from the first two terms we have 



r(i) 1 ( d2 t 3 \ 2 / 1 2 \ 

* 2 =2 [dxj ^ + pu ^ + 0^ [r U2 + pu ^) 



2 /I 



cj - 



dX2 K r ui+puiU2 

1/0/2 



+ — 1.0^ {-^dx-^ 1 

0/1/0 \\\ 1 r(0) 

"0^ \dxl U2 + dx~ 2 Ul J J - 2 0^ 2 



(107) 
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and for the third term 

d (o) _ d^ 0) dp t d^ 0) dp Ul | <9<^ 0) dpu 2 



dt dp dt dpui dt dpu 2 dt 



U 2 + u\u 2 I -r p 



d_ 

dx 2 ' 



A ,\ a fl 2 \ d n d 

+ W + 3/mi )d^ Ul + W + pu ^ ) ^ + 2puiU2 dx~ 2 Ul 

9 (fl 2 \ d fl ,\ d „ 8 



2 d 2 9 „ d 

+ P U 2 n U l + P U 1^ M 2 + 2pUiU 2 - U 2 

ox 2 oxi dx 2 



Combining all three terms we have 



(108) 



rW f"- 1 A ( 9 ( 2 9 

-sJUr Uw^ 1 

-^(-^(^ 2 + ^))) (1 ° 9) 
The final macroscopic equations for this particular lattice and quasiequilib- 
rium then are 

9 

-p = -D x pu 



dt 

2-w / <9 / 2 <9 \ d f 1 f d d 



9 9 (P, 2\ 9 



2lu \dxi \ 3 P dxi Ul ) dx 2 \ 3 P \dx\ U2 ^ dx 2 Ul 

(110) 

and again the second momentum density can be found by reflection. Once 
again many of the disippative terms vanish in the continuous velocity system. 



8 Macroscopic Stability 

In the previous sections we have demonstrated the discrete velocity systems 
studied do not recover the exact macroscopic dissipativc dynamics of the 
continuous system. We are now concerned with the stability of the discrete 
dynamics under a short wave perturbation. In each example we are concerned 
with the stability of the linear part of the dynamics (as calculated above) only. 



28 



D.J. Packwood et al 



8.1 The athermal 1-D model 

We consider perturbations by a Fourier mode around a constant flow, that is 
we write 

p = p a + Ae^ xt+KX ^ 
u = u + Be^ Xt+KX ^ 

We combine this with a composite coefficient for the first order part 

" = ^ <" 2 > 

Substituting these into the macroscopic equations and with some rearrange- 
ment for the u term we have 



AX = —poBn — uqAk 



BX = -J— An - u Bk - v f—Ain 2 + ( Zul - ^] Bin' 
We take eigenvalues of the matrix 



(113) 



— UqK —pQK 
gi^-K — V^-IK 2 —UqK — V (3Mq — |) ^ 



(114) 



which give us two values for A 
( 3 2 . 1 . 

A = K I —Mo — -VUqIK + -VIK 



±\J vu^ik - ^v 2 u\k 2 + v 2 u\k 2 - ^ 2 k 2 + ^ (115) 

In order for the manifold to remain bounded in time we investigate parame- 
ters which give 9(A) > 0. We begin by checking aymptotics of two parameters, 
for large n we have 



™ 2 J*±y- (§«8+J 

0. h'tr [ — 3uq + ^ J i 



(116) 



and for large uq 



^ f 3 2- / 9 2 4 2 \ 



2„-„2 



= 0, —3vu in 
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We can see from this that for non-zero re the first condition that should be 
satisfied for stability is Uq < 2/9, for large uq it is necessary for re to equal 0. 
Additionally, stability is absolutely contingent on the composite coefficient v 
being positive, this is the dual condition that time steps are positive and that 
relaxation parameter of the collision u> is in the interval < uj < 2 (repeated 
steps of the collision integral in isolation go towards the quasiequilibrium) . 
In the case that either v is negative or that uq is outside the given region, 
the magnitude of the Fourier perturbation will grow exponentially causing a 
rapid divergence from the constant flow. 

We can confirm these results numerically by plotting the contours of the 
two eigenvalues equal to zero. In fact in Figure© we additionally plot con- 
tours below zero to show the decay from stability. 



H 




! 

i 


0.5 

K o 

-0.5 


1 




0.5 

K 

-0.5 


H 





□ 1 2 3 -3 -2 -1 1 2 3 -3 -2 -1 

u a u a tio 



Fig. 2. The first two figures show stability for each of the two eigenvalues 
in the D1Q3 system with v = 1, the third figure plots the minimum of the 
two. Contours are plotted at 3(A) = (—0.3, —0.2, —0.1, 0), the yellow region 
and it's boundary indicate the stable region and, the other colours, the decay 
from stability. 



8.2 The athermal 2-D model 



We extend the stability analysis from the one dimensional case with a per- 
turbation in the additional space direction. The perturbed system is given 

by, 

p = po + Ae l{xt+KlXl+K2X2) 

l*i = uio + B 1 e l( - xt+Klxl+K2X2 '> (118) 
u 2 = u 20 + B 2 e l( - xt+KlXl+K2X2 1 
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In this case we investigate the short wave asymptotics as |fti|, [/ca j — > oo. The 
eigenvalues of the system under such conditions are 

f 1 3 a V 2 fl 3 a \ . 2 . 

2 ^ "VK 2 + f ^ — u 2 2 ^ ivn 2 , - 2iisu l0 u2 K 1 K 2 - 

(119) 



In the 1-D examples all terms were in even powers of k whereas in this case 
there are cross terms in the product Ki«2- Because of this it is necessary 
to consider the different permutations of signs for these terms. Since the 
condition that the third eigenvalue imposes is weaker than that of the the 
first two, which are equivalent, it is sufficient to find the region of stability 
using just one of these. Again assuming that the coefficient v is positive, the 
region is given by parameters satisfying the two conditions. 

1 3 2 \ fl 3 2 \ n 
3 ~ 2 Ul ° J + [ 3 ~ 2 U2 ° J ~ 3ui ° U2Q - 

{ ) ( ( 12 0) 

1 3 \ fl 3 \ 

3 - 2 Ul " ) + I 3 ~ 2 W2 ° ) + 3ui ° U2 ° ~ 

The plot of the region generated by these inequalities is given in Figure [SI 
Similarly to the one dimensional examaple, in the event that v is negative 
or the constant flow speed moves outside this region, the magnitude of the 
Fourier perturbation will increase exponentially in time. 

Again for specific parameters the stability can be calculated numerically. 
In the first case examine the case where K2,W2o — 0- Figure <j3j) shows the 
stability plot for the three eigenvalues and their minimum, in this case we 
see that while the eigenvalues are different from their counterparts in the 1-D 
system, the stability region is exactly the same. In Figure Q we vary k 2 and 
U20 to see what affect this has on the stability region. For a more complete 
picture we plot u± against k 2 and again plot the stability region. In Figure 
([5]) we vary k\ and M20 across the different plots. 



9 Conclusion 

In the analaysis of the continuous Boltzmann equation, the Chapman-Enskog 
procedure is known to reproduce the Navier Stokes equations [2 3 4 . This is 
achieved by a perturbation by a small parameter, the Knudsen number. At 
the zero and first orders of this parameter, repsectively, the convective and 
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Ul Ulo "10 "10 

Fig. 3. The first three figures show stability for each of the three eigenvalues 
in the athermal D2Q9 system with parameters v = 1, u-i = 0, k% = 0, the 
fourth figure plots the minimum of them. In each case the contours are plotted 
at A = (—0.3,-0.2,-0.1,0) therefore the yellow region and it's boundary 
describe the stable area. 




Uiq Uio U H\ 



Fig. 4. Stability regions for the athermal D2Q9 system. Contours are plot- 
ted at 3(A) = (—0.3,-0.2,-0.1,0), the yellow region and it's boundary 
indicate the stable region and, the other colours, the decay from stability. 
The parameters are v = 1 and additionally a) k 2 = — 0.1, U20 = — 0.5;b) 
«2 = —0.1,W2o = 0;c) k-2 = -0.1,Ti2 = 0.5; d) n 2 = 0,u 2o = -0.5; 
e) k 2 = 0,u 2o = 0; f) K 2 = 0,u 2o = 0.5; g) k 2 = 0.1, u 2Q = -0.5;h) 
k 2 = 0.1, M 2 o = 0;i) n 2 = 0.1, u 20 = 0.5. 



diffusive dynamics appear. At higher orders which were not discussed here 
the Burnett equations arise. 

The discrete Boltzmann schemes studied here are defined by the require- 
ment that the Euler equations are recovered at the zero order. In common 
with the continuous scheme, dissipative terms arise at the first order, however 




Fig. 5. Stability regions for the athermal D2Q9 system. Contours are plot- 
ted at 3(A) = (—0.3,-0.2,-0.1,0), the yellow region and it's boundary 
indicate the stable region and, the other colours, the decay from stability. 
The parameters are v = 1 and additionally a) Ki = — 0.1, U20 = — 0.5;b) 
Ki = — 0.1, M20 = 0;c) K\ = — 0.1, U20 = 0.5; d) k\ = 0,U2o = —0.5; 
e) ki = 0,u 2 o = 0; f) k x = 0,u 20 = 0.5; g) n x = 0.1, u 20 = -0.5;h) 
ki = 0.1, m 2 o = 0;i) ki = 0.1, U20 — 0.5. 




Fig. 6. Stability regions for the athermal D2Q9 system. Contours are plotted 
at 9(A) = (—0.3, —0.2, —0.1, 0), the yellow region and it's boundary indicate 
the stable region and, the other colours, the decay from stability. The param- 
eters are v = 1 and |«i|, \ — > 00. 



in the discrete case there appear additional viscous terms. In parallel with 
Goodman and Lax[6] we view the additional dissipative part of the fluid as a 
direct consequence of the discrete scheme used. In this work we have used the 
idea of invariant manifolds [7| to calculate the macroscopic dynamics arising 
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from discrete time Boltzmann schemes. This technique is based on an expan- 
sion in a different small parameter, the time step e. Dynamics at the zero 
and first orders again correspond to the conservative and dissipative parts of 
a fluid respectively. Although in this work we calculate these dynamics up 
to the first order only, the methodology can be extended to calculate higher 
order systems. 

To compute a solution to the Boltzmann system it is also necessary to dis- 
cretize velocity space. We have presented two alternative modes of thought 
to reason why this should be, which produce equivalent systems. We have 
calculated the exact macroscopic dynamics up to first order of two common 
discrete velocity schemes, and their continuous counterparts. Although the 
dynamics of these two schemes match at the zero order, in the discrete ve- 
locity case additional erroneous terms arise at the first order. Such errors 
might be expected due to the way the quasi-equilibria in the discrete case 
are defined. If we view the discrete velocity summation as a quadrature ap- 
proximation to the continuous velocity integral, then we should expect an 
error of integration. At the zero order we find no such error. This is due to an 
equilibrium being constructed specifically that the zero order moments are 
calculated exactly. This equilibrium consists of merely the first three terms 
of the Taylor expansion of the continuous equilibrium about the zero mo- 
mentum position. It should, perhaps then, be no surprise that the dissipative 
dynamics in the discrete system approach those of the continuous system 
only in the limit of momentum going to zero. 

Finally we perform a stability analysis of the linear part of the macroscopic 
dynamics of the discrete velocity schemes under a short wave perturbation. 
In common with other authors using similar Fourier techniques [17 18], and 
with our own earlier assumptions, we find that two lattice parameters are 
critical for stability. These are the time step e which must be positive, and 
the relaxation parameter w, which must be chosen for non-zero flow speed in 
the interval (— 1, 1). We also analytically and graphically give the permissible 
range of macroscopic quantities for stability. For the athermal systems study 
the density p can be any value, whereas the momentum u should be within 
an area centered around the zero point. The exact shape of this region is 
determined by the choice of velocity discretization. 
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Appendix: First order populations for the D2Q9 lattice 
with standard polynomial quasi-equilibria 
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